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Abstract.  A  probability-one  homotopy  algorithm  for  solving  nonsmooth  equations  is  described. 
This  algorithm  is  able  to  solve  problems  involving  highly  nonlinear  equations,  where  the  norm  of  the 
residual  has  non-global  local  minima.  The  algorithm  is  based  on  constructing  homotopy  mappings 
that  are  smooth  in  the  interior  of  their  domains.  The  algorithm  is  specialized  to  solve  mixed  comple¬ 
mentarity  problems  through  the  use  of  MCP  functions  and  associated  smoothers.  This  specialized 
algorithm  includes  an  option  to  ensure  that  all  iterates  remain  feasible.  Easily  satisfiable  sufficient 
conditions  are  given  to  ensure  that  the  homotopy  zero  curve  remains  feasible,  and  global  convergence 
properties  for  the  MCP  algorithm  are  developed.  Computational  results  on  the  MCPLIB  test  library 
demonstrate  the  effectiveness  of  the  algorithm. 
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1.  Introduction.  The  primary  attraction  of  homotopy  algorithms  is  that  they 
are  able  to  reliably  solve  systems  of  equations  involving  highly  nonlinear  functions, 
where  the  norm  of  the  residual  may  have  non-global  local  minima.  This  is  because, 
unlike  line  search  or  trust  region  methods,  homotopy  methods  do  not  rely  on  descent 
of  a  merit  function.  Instead,  they  work  by  following  a  path,  which  under  certain 
weak  assumptions  is  known  to  lead  to  a  solution.  Standard  probability-one  homotopy 
algorithms  require  that  the  system  of  equations  involves  only  smooth  ( C 2)  functions. 
This  paper  proposes  a  probability-one  homotopy  algorithm  for  solving  nonsmooth 
systems  of  equations,  and  specializes  this  algorithm  to  solve  mixed  complementarity 
problems.  The  algorithm  uses  smoothing  functions  to  construct  a  homotopy  mapping 
that  is  C2  in  the  interior  of  its  domain.  This  allows  the  zero  curve  of  the  homotopy 
mapping  to  be  tracked  using  software  from  the  HOMPACK90  suite  of  homotopy  codes 
[24] .  A  preliminary  version  of  this  algorithm  was  presented  at  the  Second  International 
Conference  on  Complementarity  Problems  [5].  The  algorithm  proposed  here  has  two 
significant  improvements:  first,  a  new  end  game  strategy,  which  makes  better  use 
of  available  information  about  the  behavior  of  the  homotopy  zero  curve;  second,  an 
option  for  mixed  complementarity  problems  that  ensures  that  all  iterates  generated 
by  the  algorithm  are  feasible.  This  is  important  because  many  applications  involve 
functions  that  are  not  defined  outside  of  the  feasible  region.  For  the  case  of  mixed 
complementarity  problems,  new  convergence  results  are  presented,  which  establish 
easily  satisfiable  sufficient  conditions  to  ensure  that  the  homotopy  zero  curve  always 
remains  strictly  feasible. 

In  order  to  describe  the  algorithm,  a  significant  amount  of  background  material 
is  needed.  This  is  given  in  Section  2,  which  discusses  notation,  nonsmooth  equations, 
a  generalized  Newton  method  for  nonsmooth  equations  (which  will  be  used  in  the  end 


*Department  of  Mathematics,  University  of  Colorado  at  Denver,  Denver,  CO,  802713364 
(sbillups@carbon.cudenver.edu),  research  partially  supported  through  NSF  Grant  DMS-9973321. 

'Departments  of  Computer  Science  and  Mathematics,  Virginia  Polytechnic  Institute  and  State 
University,  Blacksburg,  VA  24061-0106,  (ltw@cayuga.cs.vt.edu),  research  partially  supported  by 
AFOSR  Grant  F496320-99-1-0128  and  NSF  Grant  DMS-9625968. 


1 


Report  Documentation  Page 

Form  Approved 

OMB  No.  0704-0188 

Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 

VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  OMB  control  number. 

1.  REPORT  DATE 

2006 

2.  REPORT  TYPE 

3.  DATES  COVERED 

00-00-2006  to  00-00-2006 

4.  TITLE  AND  SUBTITLE 

5a.  CONTRACT  NUMBER 

A  Probability-One  Homotopy  Algorithm  for  Nonsmooth  Equations  and 

5b.  GRANT  NUMBER 

iviiAcu  riuuicni!) 

5c.  PROGRAM  ELEMENT  NUMBER 

6.  AUTHOR(S) 

5d.  PROJECT  NUMBER 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

University  of  Colorado  at  Denver, Center  for  Computational 

Mathematics, PO  Box  173364, Denver, CO, 80217-3364 

8.  PERFORMING  ORGANIZATION 

REPORT  NUMBER 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

10.  SPONSOR/MONITOR'S  ACRONYM(S) 

11.  SPONSOR/MONITOR'S  REPORT 
NUMBER(S) 

12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited 

13.  SUPPLEMENTARY  NOTES 

14.  ABSTRACT 

see  report 

15.  SUBJECT  TERMS 

16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION  OF 
ABSTRACT 

18.  NUMBER 
OF  PAGES 

19 

19a.  NAME  OF 
RESPONSIBLE  PERSON 

a.  REPORT 

unclassified 

b.  ABSTRACT 

unclassified 

c.  THIS  PAGE 

unclassified 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std  Z39-18 


game),  probability-one  homotopy  methods,  complementarity  problems,  and  smooth¬ 
ing  functions.  Section  3  describes  a  probability-one  homotopy  algorithm  for  nons¬ 
mooth  equations.  This  algorithm  is  then  specialized  to  solve  mixed  complementarity 
problems  in  Section  4.  Section  5  addresses  implementation  details  and  computational 
results,  and  Section  6  concludes. 

2.  Background. 

2.1.  Notation.  When  discussing  vectors  and  vector-valued  functions,  subscripts 
are  used  to  indicate  components,  whereas  superscripts  are  used  to  indicate  the  itera¬ 
tion  number  or  some  other  label.  In  contrast,  for  scalars  or  scalar-valued  functions, 
subscripts  refer  to  labels  so  that  superscripts  can  be  used  for  exponentiation.  The 
vector  of  all  ones  is  represented  by  e. 

Unless  otherwise  specified,  ||-||  denotes  the  Euclidean  norm.  For  a  set  C  C  IR", 
7 rc(^)  represents  the  orthogonal  projection  (with  respect  to  the  Euclidean  norm)  of  x 
onto  C .  The  symbol  JR+  refers  to  the  nonnegative  real  numbers.  The  extended  real 
numbers  are  denoted  by  ]R,  :=  IR, |J{ — oo,  +oo}. 

Real-valued  functions  are  denoted  with  lower-case  letters  like  /  or  (f)  whereas 
vector- valued  functions  are  represented  by  upper-case  letters  like  F  or  $.  For  a 
function  F  :  C  C  IR"  — >  IR"1,  VF( x)  is  the  mxn  matrix  whose  i,j th  element 
is  dFi(x)/dxj.  Let  D  C  IRm.  Then  F~1(D)  is  the  set-valued  inverse  defined  by 
F~l(D)  :=  {ar|  F(x)  G  Dj. 

Given  a  function  F  :  IR"  — >  1R”\  the  directional  derivative  of  F  at  x  in  the 
direction  d  is  denoted  by  F'(x\  d)  :=  limqo  (F(x  +  td )  —  F(x))/ 1,  provided  the  limit 
exists. 

2.2.  Nonsmooth  equations.  This  paper  is  concerned  with  solving  equations 
of  the  form  F(x)  =  0,  where  the  function  F  :  IR"  — ►  IR"  is  locally  Lipschitzian, 
but  not  necessarily  continuously  differentiable.  Such  nonsmooth  equations  provide  a 
unifying  framework  for  the  study  of  many  important  classes  of  problems,  including 
constrained  optimization,  finite-dimensional  variational  inequalities,  complementarity 
problems,  equilibrium  problems,  generalized  equations,  partial  differential  equations, 
and  fixed  point  problems.  The  following  definitions  will  be  used  throughout  the  paper. 

By  Rademacher’s  theorem,  since  F  is  locally  Lipschitzian,  it  is  differentiable 
almost  everywhere.  Let  Dp  be  the  set  where  F  is  differentiable.  Define  the  B- 
subdifferential  by 


The  Clarke  subdifferential  dF(x)  is  the  convex  hull  of  dBF(x). 

F  is  said  to  be  semismooth  [19]  at  x  if  it  is  directionally  differentiable  at  x  and 
for  any  V  G  dF(x  +  h),  h  — >  0, 


dBF(x )  :=  V 


3{ccfc}  — >  x,xk  G  Dp,  with 


V  =  lim  VE(xfe) 

k—>  oo 


Vh-  F'(x;h)  =o(\\h\\). 

F  is  said  to  be  strongly  semismooth  [10]  if  additionally, 

Vh  —  F'(x;  h)  =  0(\\h\\2). 

A  semismooth  function  F  :  IR"  — » IR"  is  BD-regular  at  x  if  all  elements  in  dBF{x ) 
are  nonsingular,  and  F  is  strongly  regular  at  x  if  all  elements  in  dF(x)  are  nonsingular. 
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2.3.  Newton’s  method  for  nonsmooth  equations.  One  approach  to  solving 
the  nonsmooth  equation  F(x)  =  0  is  a  generalization  of  Newton’s  method  to  semis¬ 
mooth  equations,  which  was  proposed  by  Qi  [19].  Qi’s  method  is  used  together  with 
an  Armijo  line  search  in  the  end  game  of  the  homotopy  algorithm  proposed  here.  Qi’s 
algorithm,  which  is  discussed  in  detail  in  [3],  is  shown  in  Figure  2.1.  9  in  this  algorithm 
is  the  merit  function  defined  by  6{x)  :=  ^F(x)T  F(x).  Theorem  2.1,  which  is  restated 
from  [19]  and  [10],  shows  that  this  algorithm  has  the  same  fast  local  convergence 
properties  as  the  standard  (smooth)  Newton’s  method  under  natural  generalizations 
of  the  standard  assumptions. 


Fig.  2.1.  Generalized  damped  newton  method 


Step  1  [Initialization]  Select  line  search  parameters  a,  a  £  (0,1),  a  positive 
integer  mmax,  a  starting  point  x°  £  ]R",  and  a  stopping  tolerance 
tol.  Set  k  =  0. 

Step  2  [Direction  generation]  Choose  Vk  £  dBF(xk).  If  Vk  is  singular,  stop, 
returning  the  point  xk  along  with  a  failure  message.  Otherwise  choose 
the  direction 

(2.1)  dk  =  -(Cfe)"1F(  xk). 

Step  3  [Step  length  determination]  Let  m/c  be  the  smallest  nonnegative  in¬ 
teger  m  <  rrimax  such  that 

(2.2)  9(xk  +  amdk)  -  6{xk)  <  -aam6{xk). 

If  no  such  mi.  exists,  stop;  the  algorithm  failed.  Otherwise  set  xk+1  = 
xk  +  am"dk. 

Step  4  [Termination  check]  If  9{xk+1)  <  tol  stop,  returning  the  point  xk+1. 
Otherwise,  return  to  step  2,  with  k  replaced  by  k  +  1. 


Theorem  2.1.  Suppose  that  x*  is  a  solution  of  F(x)  =  0  and  that  F  is  semis¬ 
mooth  and  BD-regidar  at  x* .  Then  the  iteration  method  defined  by  xk+1  =  xk  +  dk , 
where  dk  is  given  by  (2.1)  is  well  defined  and  convergent  to  x*  Q-superlinearly  in  a 
neighborhood  ofx*.  If  F  is  strongly  semismooth  atx* ,  the  iteration  sequence  converges 
to  x*  Q-quadratically. 

One  consequence  of  this  local  convergence  theorem  is  that  within  a  neighborhood 
of  a  BD-regular  solution  x*,  the  line  search  criterion  (2.2)  will  be  satisfied  by  = 
0.  Thus,  the  inner  algorithm  will  take  full  Newton  steps  and  achieve  the  fast  local 
convergence  rates  specified  by  the  theorem. 

The  damped  Newton  method  described  above  works  very  well  when  started  near 
a  solution,  or  when  applied  to  problems  that  are  nearly  linear  in  the  sense  that  their 
merit  functions  do  not  contain  local  minima  that  are  not  solutions. 

For  highly  nonlinear  problems,  the  damped  Newton  method  tends  to  fail  without 
a  carefully  chosen  starting  point.  The  reason,  of  course,  is  that  unless  started  close 
to  a  solution,  the  iterates  may  converge  only  to  a  local  minimum  of  the  merit  func¬ 
tion.  This  motivates  the  consideration  of  homotopy  methods,  which  are  truly  globally 
convergent. 
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2.4.  Homotopy  methods.  The  main  theory  underlying  the  present  homotopy 
method  is  summarized  in  the  following  proposition  from  [5] .  This  proposition  is  similar 
to  results  presented  in  [20]  and  [8,  Theorem  2.4];  however,  it  does  not  assume  F  itself 
to  be  differentiable.  The  path  ya  defined  in  the  proposition  “reaches  a  zero  of  F”  in 
the  sense  that  it  contains  a  sequence  {(Afc,a;fe)}  that  converges  to  (1,5:),  where  a;  is  a 
zero  of  F. 

Proposition  2.2.  Let  F  :  IR"  — »  IR"  be  a  Lipschitz  continuous  function  and 
suppose  there  is  a  C 2  map 


p:TRmx  [0, 1)  x  IR"  ->  IR” 


such  that 

1.  X7p(a,X,x)  has  rank  n  on  the  set  /?_1({0}), 

2.  the  equation  pa(t),x)  =  0,  where  pa(X,x)  :=  p(a,X,x),  has  a  unique  solution 
xa  €  IR”  for  every  fixed  a  £  IRm, 

3.  Va;pa(0,a;a)  has  rank  n  for  every  a  €  lRm, 

4-  p  is  continuously  extendible  (in  the  sense  of  Buck  [6])  to  the  domain 
IRm  x  [0, 1]  x  ]R”?  and  pa(  1,  x)  =  F(x)  for  all  x  £  IR"  and  a  £  JRm,  and 
5.  7 a,  the  connected  component  o/p“1({ 0})  containing  (0,a:a),  is  bounded  for 
almost  every  a  £  lRm. 

Then  for  almost  every  a  £  IRm  there  is  a  zero  curve  ya  of  pa,  along  which  Vpa  has 
rank  n,  emanating  from  (0,xa)  and  reaching  a  zero  x  of  F  at  A  =  1.  Further,  ja  does 
not  intersect  itself  and  is  disjoint  from  any  other  zeros  of  pa  ■  Also,  if  ”fa  reaches  a 
point  (1,5:)  and  F  is  strongly  regidar  at  x,  then  ya  has  finite  arc  length. 

Because  y0  is  a  smooth  curve,  it  can  be  parameterized  by  its  arc  length  away 
from  (0,a:a).  This  yields  a  function  (A(s),x(s)),  representing  the  point  on  7a  of  arc 
length  s  away  from  (0,a:a). 

The  construction  of  a  globally  convergent  probability-one  homotopy  algorithm 
entails:  (1)  constructing  a  map  p  according  to  Proposition  2.2,  (2)  choosing  a  £  IRm, 
(3)  finding  xa  solving  pa(0,x)  =  0,  and  (4)  tracking  y0  starting  from  (0,a:a)  until 
A  =  1.  Assuming  an  appropriate  p  exists,  the  theory  guarantees  that  for  almost  all  a 
(in  the  sense  of  Lebesgue  measure),  y0  exists  and  leads  to  a  solution,  hence  the  term 
“probability-one” . 

A  simple  (and  occasionally  useful  in  practice)  homotopy  mapping  is  p  :  IR”  x 
[0, 1)  x  IR"  ->  IR”  given  by 

(2.3)  p(a,  A,  x )  :=  A F(x)  +  (1  —  A)(x  —  a). 

If  F  is  C2  then  p  trivially  satisfies  properties  (1),  (2),  (3),  and  (4)  but  not  necessarily 
(5)  of  Proposition  2.2.  The  following  theorem  gives  conditions  on  F  under  which  the 
fifth  condition  is  satisfied.  This  result  will  be  generalized  to  nonsmooth  functions  in 
Theorem  3.2. 

Theorem  2.3.  [22]  Let  F  :  IR"  — >  IR"  be  a  C2  function  such  that  for  some 

x  £  IR"  and  r  >  0, 

(2.4)  (a:  —  x)T F(x)  >  0  whenever  \\x  —  i||  =  r. 

Then  F  has  a  zero  in  a  closed  ball  of  radius  r  about  x,  and  for  almost  every  a  in  the 
interior  of  this  ball  there  is  a  zero  curve  y0  of 

pa( A,x)  :=  A F(x)  +  (1  -  X)(x  -  a), 
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along  which  \7pa(\,x)  has  full  rank,  emanating  from  (0,a)  and  reaching  a  zero  x  of 
F  at  A  =  1.  Further,  7a  has  finite  arc  length  if\7F(x)  is  nonsingular. 

The  actual  statement  of  the  theorem  in  [22]  fixes  x  =  0.  However,  the  proof 
can  be  modified  trivially  to  yield  the  more  general  theorem  above.  (See  the  proof  of 
[5,  Theorem  2.11]  for  the  necessary  modifications).  It  is  interesting  to  note  that  in 
many  applications,  (2.4)  holds  for  all  r  sufficiently  large  (not  just  for  some  fixed  r). 
This  makes  the  choice  of  x  irrelevant.  Furthermore,  in  such  cases,  a  can  be  chosen 
arbitrarily,  (instead  of  from  some  neighborhood  of  x ) ,  thus  making  the  method  truly 
globally  convergent  (with  probability  one). 

(2.4)  will  be  referred  to  as  the  global  monotonicity  property.  If  a  G2  function 
F  possesses  this  property,  these  theoretical  results  have  some  profound  implications: 
the  guaranteed  existence  of  a  path  between  almost  any  starting  point  and  a  solution 
x  to  F(x)  =  0,  which  has  finite  arc  length  if  rankV-F(ir)  =  n.  In  theory,  to  find  a 
solution,  one  must  simply  follow  the  path  to  a  point  of  ya  where  A  =  1.  In  practice, 
however,  the  task  of  constructing  a  p  for  which  ya  is  short  and  smooth  is  very  difficult, 
although  this  has  been  done  for  large  classes  of  problems. 

Several  packages  exist  to  solve  root  finding  problems  using  homotopy  techniques 
[24] .  The  implementation  here  uses  the  routine  STEPNX  from  the  HOMPACK90  suite 
of  software  [23]  [24,  Section  3],  which  tracks  the  zero  curve  of  a  homotopy  mapping 
specified  by  the  user. 

2.5.  Complementarity  problems.  Given  a  continuously  differentiable  func¬ 
tion  G  :  IR"  — >  IR",  the  nonlinear  complementarity  problem  NCP(G)  is  to  find  some 
x  €  R"  so  that 

(2.5)  0  <  x  L  G{x)  >  0, 

where  x  T  G(x)  means  that  xTG(x )  =  0. 

Given  a  rectangular  region  JBj.u  :=  C  IR,"  defined  by  two  vectors, 

l  and  u  in  IR.”  where  —oo<l<u<oo,  and  a  function  G  :  IR”  — ►  IR”,  the 
mixed  complementarity  problem  MCP(G,lBz,u)  is  to  find  an  x  G  IB i,u  such  that  for 
each  i  G  {l,...,n},  either  1)  Xi  =  li  and  Gi(x )  >  0,  2)  Gifx )  =  0,  or  3)  Xi  = 
Ui  and  Gi(x )  <  0.  This  is  equivalent  to  the  condition  that  mid (x  —  l,x  —  u,G(x))  =  0, 
where  mid  represents  the  componentwise  median  function.  When  these  conditions 
are  satisfied,  write  G(x)  A.  x  and  say  that  x  is  complementary  to  G(x).  Assume 
henceforth  that  G  is  G2. 

It  is  well  known  that  NCP(G)  can  be  reformulated  as  a  system  of  equations.  This 
was  first  shown  by  Mangasarian  [16].  An  excellent  review  of  reformulations  of  NCP 
can  be  found  in  [18].  To  discuss  such  reformulations  requires  several  definitions,  which 
are  equivalent  to  the  NCP  function  and  the  BVIP  function  defined  in  [18]: 

Definition  2.4.  A  function  (j)  :  ]R2  — »  IR,  is  called  an  NCP  function  provided 
<p(a ,  b)  =  0  if  and  only  if  min(a,  b)  =  0. 

Definition  2.5.  A  function  ip  :  IR (J{  —  oo}  x  1R,U{°°}  x  IR2  — ►  IR,  is  called  an 
MCP  function  provided  ip (l,  u,  a,b)  =0  if  and  only  z/mid(a  —  l,a  —  u,b)  =  0. 

It  is  useful  to  further  distinguish  NCP  and  MCP  functions  according  to  their 
orientations: 

Definition  2.6.  An  NCP  function  <p  is  called  positively  oriented  if  for  all  a,b  G 
IR, 


sign(^(a,  b))  =  sign(min(a,  b )). 
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An  MCP  function  ip  is  called  positively  oriented  if 


sign  {ip  {l,  u,  a,b))  =  sign(mid(a  —  l,  a  —  u,  b)) 

for  all  l  €  ]R  U{— °o},  u  G  ]R,  U{oo},  l  <  u,  and  a,b  €  ]R. 

An  NCP  function  that  has  been  very  popular  recently  is  the  Fischer-Burmeister 
function  [13]  (p  :  IR2  — » IR,  defined  by 

(2.6)  < f>FB(a ,  b )  :=  a  +  b  —  \/ a2  +  b2. 

It  is  easily  seen  that  <pFB(a ,  b)  =  0  if  and  only  if  0  <  a  _L  b  >  0.  Thus,  by  defining  the 
function  F  :  IR"  — ■>  IR”  by 

(2.7)  Fi(x)  :=  <pFB(xi,Gi(x)), 

it  is  clear  that  x  €  IR"  solves  NCP(G)  if  and  only  if  F(x)  =  0. 

While  (pFB  is  not  differentiable  at  the  origin,  (<pFB)2  is  continuously  differentiable 
everywhere.  This  property,  together  with  the  fact  that  <pFB  is  semismooth,  makes  this 
reformulation  well  suited  for  use  in  globalization  strategies  for  nonsmooth  Newton- 
based  methods  (see,  for  example,  [9]). 

Given  a  positively  oriented  NCP  function  cp,  and  the  convention  that  </>(oo,6)  = 
linia^oo  <p(a,  b)  and  </>(a,  oo)  =  lim^oo  4>{a ,  b),  an  MCP  function  ip  can  be  constructed 
using  the  following  formula,  first  proposed  in  [1]: 


(2.8)  ip  (l,  u,  a,  b)  :=  cp(a  —  l,  —<p(u  —  a,  —b)). 

Constructing  the  function  F  :  JR"  — >  R"  by 

(2.9)  Fi (xj  . —  1p(li ,  Ui ,  Xi ,  Gi (x)) , 

yields  a  reformulation  of  the  MCP(G,  IB;  „);  F( x)  =  0  if  and  only  if  x  is  a  solution  to 
MCP(G,  JB;jtl)  [2], 

Note  that  for  the  Fischer-Burmeister  function,  lim^oo  cpFB(a,  b)  =  b  and 
lim^oo  (pFB{a,  b)  =  a.  Thus,  for  the  MCP  case,  if  7;  is  finite,  and  Ui  =  oo,  then 
Fi(x)  =  cpFB{xi  —  li,Gi{x ));  if  is  finite  and  li  =  — oo,  then  F.^x)  =  —<pFB(ui  — 
Xi,—Gi(x));  and  if  neither  bound  is  finite,  Fi(x)  =  Gi(x). 

2.6.  Smoothing  operators.  Consider  the  system  F(x)  =  0  where  F  is  a  nons¬ 
mooth  function,  and  suppose  there  exists  a  family  of  functions  F M  parameterized  by  a 

smoothing  parameter  p  so  that  lim  F M  =  F  in  some  sense.  Under  suitable  conditions, 

MO 

the  solutions  to  the  systems  F^(x)  =  0  converge  to  a  solution  to  F(x)  =  0  along  a 
smooth  trajectory  [7]. 

Definition  2.7.  Given  a  nonsmooth  continuous  function  <p  :  IRP  — >  JR,  a 
smoother  for  <p  is  a  continuous  function  <p  :  ]RP  x  ]R+  — » IR  such  that 

1.  <p{x,  0)  =  <p(x),  and 

2.  (p  is  continuously  differentiable  on  the  set  ]RP  x  1R.++. 

If  <j>  is  G2  on  ]RP  x  1R.++,  call  <j>  a  G2-smoother. 

For  convenience,  define  (pffx)  :=  <p(-,p).  To  define  smoothers  for  functions  F  : 
IR”  — » IR”,  say  that  F M  :  IR”  x  1R.+  — >  IR"  is  a  smoother  for  F  if  for  each  i  G  {1 . . .  n}, 
Ff  is  a  smoother  for  Ft. 
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In  the  case  of  complementarity  problems,  the  NCP  functions  and  MCP  functions 
generally  have  well  understood  nonsmoothness  structure,  so  C2-smoothers  for  these 
functions  can  usually  be  easily  constructed.  As  an  example,  the  following  C^-smoother 
for  the  Fischer-Burmeister  function  was  proposed  by  Kanzow  [15]: 

(2.10)  b,fi)  :=  a  +  b  —  y/a2  +  b2  +  2/z, 

The  following  smoother  is  more  useful  here,  since  its  partial  derivative  with  respect 
to  [i  is  bounded  near  the  origin. 

(2.11)  (j>BW (a,  b,n)  :=  a  +  b  -  a2  +  b2  +  /x2, 

Given  a  smoother  <j>  for  a  NCP  function  <f>  and  the  convention  that  <j>(oo,b,  n)  = 
linia^oo  b,  /j.)  and  <p(a,oo,  fi)  =  lim^oo  b,  fi),  a  smoother  %p  for  the  MCP 
function  ip  defined  by  (2.8)  can  be  constructed  according  to  the  formula: 

(2.12)  r/>(Z,  u,  a,  b ,  n)  :=  </>p(a  -  l,  -0M(u  -  a,  - b )). 

Smoothers  for  (2.7)  and  (2.9)  are  then  given,  respectively,  by 


(2.13)  F?{x)  :=  (j>^,{xi,Gi{x)),  and 

(2.14)  Fi(x)  '■=  i’n(k,ui,Xi,Gi(x)). 

Note  that  for  the  smoother  defined  by  (2.11),  lim^oo  <f)BW (a,b,  n)  =  6,  and 
lim^oo  (j)BW(a,b,  n)  =  a.  Thus,  for  the  MCP  case,  if  ul  —  oo  and  lj  is  finite,  then 
F?{x)  =  (pBW(xi  —  U,  Gi(x),  /z);  if  Ui  is  finite  and  lt  =  — oo,  then  F^{x)  =  —(j>BW{ui  — 
Xi,  —  Gi(x),(j,);  and  if  neither  bound  is  finite,  then  F^x)  =  Gi(x). 

3.  The  algorithm.  This  section  summarizes  the  probability-one  homotopy  al¬ 
gorithm  for  solving  nonsmooth  equations.  It  contrasts  with  an  earlier  hybrid  Newton- 
homotopy  method  described  in  [2] .  The  earlier  method  begins  by  using  a  nonsmooth 
version  of  a  clamped-Newton’s  method  to  solve  the  root  finding  problem  F{x)  =  0. 
If  the  Newton  algorithm  stalls,  a  standard  homotopy  method  is  invoked  to  solve  a 
particular  smoothed  version  of  the  original  problem,  F"M( x)  =  0,  where  /i  is  fixed.  The 
smoothing  parameter  /z  is  chosen  based  on  the  level  of  a  merit  function  on  F  at  the 
last  point  x  generated  by  the  Newton  method.  Starting  from  x,  a  homotopy  method 
is  carried  out  until  it  produces  a  point  that  yields  a  better  merit  value  than  the  pre¬ 
vious  Newton  iterate.  The  Newton  method  is  then  started  again  and  the  process 
repeats  until  a  point  is  produced  that  is  close  enough  to  a  solution  or  the  homotopy 
method  fails.  One  key  feature  of  that  hybrid  method  is  that  each  time  the  Newton 
method  stalls,  a  different  homotopy  map  is  constructed.  The  smoothing  parameter  /.t 
is  chosen  based  on  the  level  of  the  merit  function  when  the  Newton  method  stalls,  so 
the  homotopy  that  is  then  used  is 

Pa{\  x)  :=  XF^(x)  +  (1  -  A)(x  -  a). 

An  alternative  approach,  described  here,  is  to  adopt  a  pure  probability-one  ho¬ 
motopy  algorithm  by  fixing  the  homotopy  map  and  tracking  a  single  homotopy  zero 
curve  into  the  Newton  domain  of  convergence  around  a  solution.  Essentially,  the  idea 
is  to  use  a  standard  probability-one  homotopy  algorithm,  but  with  a  specially  de¬ 
signed  “end  game”  near  a  solution.  The  key  to  this  approach  is  to  define  a  homotopy 
mapping  that  couples  the  smoothing  parameter  with  the  homotopy  parameter. 
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3.1.  The  homotopy  map.  Given  a  function  F  and  an  associated  C-smoother 
F11,  construct  a  homotopy  mapping  with  where  the  smoothing  parameter  p  is 
a  function  of  the  homotopy  parameter  A  so  that  p  j  0  as  A  f  1.  If  this  homotopy 
satisfies  the  conditions  in  Proposition  2.2,  a  well  behaved  path  exists  from  almost  any 
starting  point  to  a  solution,  and  standard  curve  tracking  techniques  can  reliably  solve 
the  equation  F(x)  =0. 

Throughout  this  section,  assume  that  F  is  a  Lipschitz  continuous  function  on  ]R” 
and  that  is  a  C2-smoother  for  F.  Take  p  :  [0, 1]  — >  ]R+  to  be  a  decreasing  C2 
function  such  that  p{ A)  >  0  for  A  <  1  and  /x(l)  =  0.  For  example, 

(3.1)  p(X)  :=  a(l  —  A) 

for  some  parameter  a  >  0.  Define  the  homotopy  map  pa  :  [0, 1)  x]Rra  — » IR",  nonlinear 
in  A,  by 

(3.2)  pa{ X,x)  :=XF^x\x)  +  (l-X)(x-a) 

and  let  ya  be  the  connected  component  of  the  set  p~  1({0})  that  contains  (0,  a).  Notice 
that  this  mapping  is  a  generalization  of  (2.3),  since  if  F  is  C2,  then  :=  F  suffices. 

In  order  to  ensure  that  a  well  behaved  zero  curve  exists,  conditions  on  F  and  its 
smoother  are  required  so  that  Proposition  2.2  can  be  invoked.  The  following  weak 
assumption  on  the  smoother  will  be  useful  in  the  theory  that  follows. 

Assumption  3.1.  There  is  a  nondecreasing  function  p  :  IR.+  — >  IR+  satisfying 
lim  p(v)  =  0  such  that  for  all  x  in  IR"  and  all  v  in  JR_|_ 

i'|0 

\\Fu(x)  —  F(a;)||oo  <  r){v). 

Note  (by  [2,  Proposition  2.14])  that  if  Fv  is  constructed  by  (2.14),  with  defined 
either  by  (2.10)  or  (2.11),  Assumption  3.1  is  satisfied  with  r\(v)  :=  3\/2n  or  p{v)  :=  3v, 
respectively. 

The  following  theorem  [5,  Theorem  2.11]  is  a  generalization  of  Theorem  2.3. 
Theorem  3.2.  Let  F  :  ]Rra  — »  IR™  he  a  Lipschitz  continuous  function  such  that 
for  some  fixed  r  >  0  and  x  £  ]R” , 

( x  —  x)T F(x)  >  0  whenever  ||ar  —  x||  =  r, 

and  let  be  a  smoother  for  F  satisfying  Assumption  3.1.  Further,  suppose  that  the 
smoothing  parameter  p(X)  is  such  that 

(3.3)  r](p(X))  <  for  0  <  A  <  1 

A 

for  some  M  £  (0,r).  Then  is  bounded  for  almost  every  a  £  IR"  such  that  ||  a  —  x||  < 
r  :=  r —  M . 

A  direct  application  of  Proposition  2.2  gives  the  main  convergence  theorem. 
Theorem  3.3.  Under  the  assumptions  of  Theorem  3.2,  F  has  a  zero  in  a  closed 
ball  of  radius  r  about  x,  and  for  almost  every  a  in  the  interior  of  a  ball  of  radius  r 
about  x,  there  is  a  zero  curve  ya  of 

p(a,  A,  x)  :=  pa{ X,x)  :=  A  FM(A)(a;)  +  (1  -  A)(x  -  a), 

along  which  Vpa( X,x)  has  full  rank,  emanating  from  (0,a)  and  reaching  a  zero  x  of 
F  at  A  =  1.  Further,  •ja  has  finite  arc  length  if  F  is  strongly  regular  at  x. 

Observe  that  in  applications,  the  r  in  Theorem  3.2  can  be  arbitrarily  large,  hence 
so  can  r  =  r  —  M,  and  thus  ||a  —  x|  <  f  is  really  no  restriction  at  all. 


3.2.  Tracking  the  zero  curve.  As  discussed  in  Section  2.4,  the  zero  curve  can, 
with  probability  one,  be  parameterized  by  arc  length:  let  (A(s),x(s))  be  the  point 
on  7 a  of  arc  length  s  away  from  (0,a;a).  Tracking  the  zero  curve  involves  generating 
a  sequence  of  points  {yk}  C  JR”+1,  with  y°  =  (0,  x“)  that  lie  approximately  on  the 
curve  in  order  of  increasing  arc  length.  That  is,  yk  «  (A (sk),  x(sk)),  where  {s*,}  is 
some  increasing  sequence  of  arc  lengths. 

The  subroutine  STEPNX  from  HOMPACK90  [24]  is  used  to  handle  the  curve 
tracking.  At  each  iteration,  STEPNX  uses  a  predictor-corrector  algorithm  to  gener¬ 
ate  the  next  point  on  the  curve.  The  prediction  phase  requires  for  each  iterate  yk 
the  corresponding  unit  tangent  vector  to  the  curve,  ( y')k  «  (X'(sk),  x'(sk))-  This 
is  accomplished  by  finding  an  element  y  of  the  null  space  of  Vpa(j/fc),  and  setting 
(y')k  :=  ±77/  ||»7||,  where  the  sign  is  chosen  so  that  ( y')k  makes  an  acute  angle  with 
(y,)fc_1)  f°r  k  >  0.  On  the  first  iterate,  the  sign  is  chosen  so  that  the  first  component 
(corresponding  to  A)  of  (7/)0  is  positive. 

At  each  iteration  after  the  first,  STEPNX  approximates  the  zero  curve  with  a 
Hermite  cubic  polynomial  ck(s ),  which  is  constructed  using  the  last  two  points  yk~x 
and  yk,  along  with  the  associated  unit  tangent  vectors  ( y ')k~1  and  ( y')k .  A  step 
of  length  h  along  this  cubic  yields  the  predicted  point  wk,°  :=  c(sk  +  h).  The  first 
iteration  uses  a  linear  predictor  instead,  which  is  constructed  using  the  starting  point 
y°  and  its  associated  unit  tangent  vector. 

Once  the  predicted  point  is  calculated,  a  normal  flow  corrector  algorithm  [24]  is 
used  to  return  to  the  zero  curve.  Starting  with  the  initial  point  wk’°,  the  corrector 
iterates  wk’\j  =  1,...  are  calculated  via  the  formula  wk,i+1  :=  wk ,J'  +  j  = 
0, 1, . . .,  where  the  step  zk^  is  the  unique  minimum-norm  solution  to  the  equation 

(3.4)  Vpaiw^Z^  =  -Pa(wk'j). 

The  corrector  algorithm  terminates  when  one  of  the  following  conditions  is  satisfied: 
the  normalized  correction  step  zk':l  /(l  +  |uAJ  ||)  is  sufficiently  small;  some  maximum 
number  of  iterations  (usually  4)  is  exceeded;  or  a  rank-deficient  Jacobian  matrix  is 
encountered  in  (3.4).  In  the  first  case,  set  yk+l  :=  ur,J’,  calculate  an  optimal  step 
size  h  for  the  next  iteration,  and  proceed  to  the  next  prediction  step.  In  the  second 
case,  discard  the  point  and  return  to  the  prediction  phase  using  a  smaller  step  size  if 
possible;  otherwise,  terminate  curve  tracking  with  an  error  return.  In  the  third  case, 
terminate  the  curve  tracking,  since  rankVpa  <  n  should  theoretically  not  happen, 
and  indicates  serious  difficulty.  The  step  size  in  h  is  also  never  reduced  beyond  relative 
machine  precision. 

3.2.1.  Step  size  control.  At  each  iteration,  STEPNX  estimates  an  “optimal” 
step  size  to  be  used  in  computing  the  predicted  point.  This  calculation  is  governed  by 
several  user-defined  parameters.  Successful  termination  of  the  corrector  phase  occurs 
when  the  norm  of  the  residual  1 1  p(xt?fc’J )  1 1  is  sufficiently  small.  In  some  cases,  this  can 
happen  even  when  the  converged  point  is  not  close  to  the  true  zero  curve.  As  the 
tracking  progresses,  the  computed  points  may  slowly  drift  farther  and  farther  from 
the  zero  curve,  while  continuing  to  meet  the  criterion  on  the  norm  of  the  residual. 
Eventually,  the  iterates  may  leave  the  Newton  domain  of  attraction,  and  the  corrector 
phase  may  fail  to  converge,  no  matter  how  small  the  predictor  step  is.  To  avoid  such 
difficulties,  STEPNX  calculates  several  quantities  that  measure  the  “quality”  of  the 
step. 

The  first  quantity  is  the  contraction  factor 

||  ^fc,l ||  /||2M|| 
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which  measures  how  much  the  Newton  step  shrinks  from  the  first  corrector  iteration 
to  the  second.  The  second  quantity  is  the  residual  factor 

||Pa(wfc’1)||  /||pa(wfc’°)||  . 

The  third  quantity  is  the  distance  factor 

Ww*’1  -  yk+1\\ /\\wk’°  -  yk+1\\  , 

which  approximates  how  much  the  distance  from  the  zero  curve  shrinks  from  the 
first  iteration  to  the  second.  Since  Newton’s  method  has  quadratic  local  convergence, 
each  of  these  quantities  should  be  small  when  the  predicted  point  is  close  to  the 
zero  curve.  Through  the  use  of  input  parameters,  the  user  is  able  to  specify  ideal 
values  (lideal,  rideal,  dideal,  respectively)  for  each  of  these  quantities.  If  the 
quantities  are  smaller  than  the  ideal,  the  step  size  will  be  increased;  if  the  quantities 
are  larger  than  ideal,  the  step  size  will  be  decreased.  The  amount  of  increase  or 
decrease  is  also  controlled  by  user-defined  parameters.  Generally,  default  values  for 
all  of  these  parameters  work  very  well.  However,  occasionally,  it  is  necessary  to  choose 
more  conservative  parameter  values  in  order  to  avoid  losing  the  zero  curve. 

As  a  final  consideration,  the  default  limit  on  the  number  of  Newton  iterations  in 
the  corrector  phase  is  4  (a  HOMPACK90  parameter).  In  some  cases,  increasing  this 
limit  to  6  or  8  improved  performance. 

3.3.  The  end  game.  The  standard  homotopy  method  used  by  HOMPACK90 
concludes  the  curve  tracking  with  an  end  game  strategy  that  zeros  in  on  a  point  (A,  x) 
on  the  zero  curve  with  A  =  1.  This  end  game  strategy,  which  is  a  robust  blend  of 
secant  iterations  with  Newton  corrections,  is  begun  when  a  point  (A,  x)  is  found  on 
the  zero  curve  with  A  >  1.  However,  this  approach  requires  that  p( A,  x)  be  defined  for 
A  >  1 — a  requirement  that  is  not  desirable  here  since  the  smoother  may  not  be 

defined  for  A  >  1.  Therefore  the  standard  end  game  is  replaced  with  the  generalized 
Newton  method  given  in  Figure  2.1,  which  is  begun  while  A  <  1  still. 

The  Newton  end  game  is  invoked  when  one  of  the  following  criteria  is  satisfied: 

1.  The  point  generated  by  the  cubic  predictor  (with  step  length  h)  has  A  >  1. 

2.  A  linear  predictor  with  the  same  step  length  has  A  >  1. 

3.  The  corrector  phase  of  the  algorithm  generates  a  point  with  A  >  1. 

In  all  cases,  a  starting  point  for  the  Newton  end  game  is  the  prediction  of  where  the 
zero  curve  crosses  the  hyperplane  A  =  1.  The  precise  details  follow. 

1.  First,  try  to  find  a  point  (Ac,  xc)  for  which  the  cubic  approximation  has  Ac  =  1. 
If  this  point  occurs  within  a  step  length  shorter  than  2 h,  then  xc  will  be  the 
starting  point. 

2.  Otherwise,  find  a  point  (A;,  xl )  for  which  the  linear  approximation  has  X1  =  1. 
Then  xl  will  be  the  starting  point. 

If  the  curve  tracking  fails  for  any  reason  before  the  end  game  criteria  are  met, 
then  attempt  the  nonsmooth  Newton’s  method  with  the  starting  point  x,  where  (A,  x) 
is  the  last  point  found  on  the  zero  curve. 

The  starting  point  generated  by  the  above  procedure  is  usually  quite  good.  How¬ 
ever,  in  some  cases,  the  Newton  end  game  may  fail  to  converge.  In  that  event,  simply 
return  to  tracking  the  zero  curve,  picking  up  from  the  last  point  yk  on  7 a,  but  with 
the  step  size  (computed  by  STEPNX)  cut  in  half,  and  with  the  STEPNX  tracking 
tolerances  abserr  and  relerr  also  reduced. 

Note  that  this  approach  differs  from  the  end  game  strategy  described  in  [5],  which 
simply  invoked  the  Newton  end  game  with  a  starting  point  x  whenever  a  point  (A,  x) 
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was  found  on  the  zero  curve  with  A  sufficiently  close  to  1.  The  new  end  game  strategy 
has  two  main  advantages  over  this  earlier  approach.  First,  using  the  cubic  predictor 
to  estimate  where  the  zero  curve  crosses  A  =  1  results  is  a  significantly  more  accurate 
approximation  for  the  solution  as  a  starting  point  for  Newton’s  method.  Second, 
the  new  method  takes  better  advantage  of  available  information  in  determining  when 
to  enter  the  end  game.  Specifically,  on  difficult  problems,  the  Newton  domain  of 
convergence  near  the  final  solution  will  be  small,  so  it  is  desirable  to  track  the  zero 
curve  very  close  to  A  =  1  before  trying  Newton’s  method.  This  is  exactly  what 
happens  since,  in  this  case,  the  step  size  will  likely  be  very  small.  In  contrast,  for 
easier  problems,  larger  step  sizes  will  be  used,  and  the  end  game  will  be  started 
earlier.  Again  this  is  acceptable  because  the  Newton  domain  of  convergence  around 
the  solution  will  likely  be  large. 

In  order  to  solve  the  system  F( x)  =  0,  the  nonsmooth  Newton’s  method  requires 
that  F  be  semismooth.  If,  in  addition,  F  is  BD-regular  at  a  solution  x*,  Newton’s 
method  will  converge  superlinearly  in  some  neighborhood  about  x* .  Theoretically,  to 
use  the  homotopy  approach  and  guarantee  the  end  game’s  success,  F  should  satisfy 
the  global  monotonicity  property  and  be  strongly  regular  at  every  solution.  This 
guarantees  that  the  homotopy’s  zero  curve  crosses  the  hyperplane  A  =  1  transversally 
rather  than  tangentially,  and  ensures  that  the  zero  curve  will  have  finite  arc  length. 
For  most  homotopies  used  in  practice  in  other  contexts,  even  if  the  zero  curve  is 
tangent  to  the  hyperplane  A  =  1,  a  point  with  A  >  1  near  p~  1({0})  will  be  generated, 
and  the  usual  end  game  provided  in  HOMPACK90  will  succeed  (to  modest  accuracy, 
since  X7F(x)  is  singular). 

4.  Solving  mixed  complementarity  problems.  This  section  specializes  the 
algorithm  described  above  in  order  to  solve  mixed  complementarity  problems.  The 
approach  taken  here  is  to  reformulate  the  MCP  by  defining  the  function  F  :  JRn  — » ]Rra 
according  to  (2.8)  and  (2.9),  where  (f>  is  a  positively  oriented  NCP  function,  and 
defining  a  smoother  for  F  according  to  (2.12)  and  (2.14),  where  (j)^  is  a  smoother 
for  cj).  Once  these  functions  are  defined,  the  homotopy  algorithm  described  in  the 
previous  section  can  be  used  to  find  a  zero  of  F,  which  corresponds  to  a  solution 
of  MCP.  Because  of  the  special  structure  of  these  functions,  stronger  convergence 
results  are  possible  than  for  the  general  nonsmooth  equations  problem.  The  first 
results  presented  in  this  section  are  tailored  to  particular  choices  of  <j>  and  </>M,  namely 
the  Fischer-Burmeister  NCP  function  (2.6),  and  the  smoother  (2.11).  More  general 
results  are  given  in  Theorem  4.3  and  Corollary  4.4.  In  describing  these  results  it  will 
be  useful  to  refer  to  the  following  index  set: 

h,u  =  {*  I  -oo  <k  <Ui  <  oo  }  . 

That  is,  Ipu  is  the  set  of  indices  for  which  both  the  lower  and  upper  bounds  are  finite. 

Theorem  4.1.  Let  (j)  be  the  positively  oriented  NCP  function  in  (2.6),  and  let  <p 
be  the  smoother  for  f>  in  (2.11).  Let  if  be  defined  by  (2.8)  with  associated  smoother  if 
defined  by  (2.12).  Choose  a  £  int  IB;,u-  Let  F ^  be  defined  by  (2.14),  where  p  :  [0, 1]  — > 
]R+  is  a  decreasing  C2  function  satisfying  p(l)  =  0  and 

(4.1)  p(X)2  <  2}  \ Ui  -  af){ui  -  If)  for  all  i  £  Ipu,  A  £  (0, 1]. 

A 

Define  pa  :  [0,1)  x  ]Rn  — >  ]R"  by  (3.2),  and  let  ya  be  the  connected  component  of 
p~f 1  ({0})  containing  (0,a).  Then  y0  is  contained  in  [0,1)  x  (intJBz,«). 
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Proof.  Let  (A,  x)  be  an  arbitrary  point  on  ja.  If  A  =  0,  then  x  =  a  £  intIBz,M;  so 
assume  0  <  A  <  1.  First  suppose  that  X*  <  U  for  some  i.  Then 

0  =  pi( A,x)  =  A Ff(-X\x)  +  (1  -  A )(xi  -  a,i) 
or 

(4.2)  F?{X\x)  =  -^(x*  -  Oi)  >  0, 

A 

where  the  last  inequality  follows  from  x,;  <  li  <  at,  since  a  is  interior  to  Also 

F^X\x)  =  4>(xi  -  k,C ,M),  where  C  :=  -</>(«*  -  xi;  -G;(x),/i).  Thus, 

F^~X\x)  =  4>(xi  -  k,C,p)  <  </>(£»  -  h,  C)  <  0, 

contradicting  (4.2).  It  follows  that  every  point  (A,  x)  on  7a  satisfies  l  <  x. 

Now  suppose  Xi  >  u,  for  some  Note  that  this  implies  that  Ui  is  finite.  In  this 
case  (analogous  to  (4.2)), 

(4.3)  F^X\x)  =  -1  ~^{xi  -  at)  <  -1  ~  ^(ui  -  ai) 

A  A 

and  £  =  —(j)(ui  —  Xi,  —  Gi(x),  n)  >  0  since  /i( A)  >  0  for  A  <  1.  If  li  =  — oo,  then 

Ffl^x\x)  =  C  >  0,  contradicting  (4.3).  If  k  is  finite,  then  from  (6)  and  (11),  for  any 

ex.,  (3  £  ]R, 

(4.4)  (j>(a,  0)  >  -  ^  . 

2\/cr  + 

Then,  using  (  >  0,  Xj  >  Uj,  the  monotonicity  of  </>,  (4.4)  and  (4.1)  gives 


contradicting  (4.3).  Therefore  every  point  (A,x)  on  satisfies  l  <  x  <  u.  □ 

Note  that  if  Ii,u  is  empty,  then  the  condition  on  p(X)  in  the  above  theorem  is 
achieved  by  any  decreasing  C2  function  satisfying  p(l)  =  0.  If  J;)U  is  not  empty,  the 
condition  is  easily  achieved  by  choosing  a  deep  in  the  interior  of  the  feasible  region 
H3;,u-  For  example,  if  Ui  —  ai>  )2  (ut  —  If)  for  all  i  £  Ii  u,  then 


min  ( Ui 

ieli.v. 


k)  (1-A) 
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suffices,  since  for  0  <  A  <  1 


MA)2 


min  (ui 

ieii,u 


1 2 

If)  (1  -  A)2 


<  2 

min  (ui 

li) 

(1-A); 

i£ll,u 

<  2 

min  (uj 

-  ai){ui  -  U) 

(l-A) 

A 

The  above  theorem  has  two  important  consequences.  First,  because  7a  always 
stays  in  the  feasible  region,  it  is  possible  to  implement  the  algorithm  without  ever 
having  to  evaluate  functions  outside  of  the  feasible  region.  This  is  important  because 
many  applications  involve  functions  that  are  not  defined  outside  the  feasible  region. 
The  second  consequence  is  the  guarantee  that  when  all  bounds  are  finite,  the  zero 
curve  7a  is  bounded.  The  implications  of  this  are  stated  in  the  following  corollary. 

Corollary  4.2.  Let  cf>  and  be  defined  by  (2.6)  and  (2.11),  respectively.  As¬ 
sume  that  all  the  bounds  of  the  MCP  are  finite,  choose  k  £  (0,  y/2 )  and  take 


(4.5) 


A4  (A) 


K 


min(wi 

.  i 


h)  (l-A), 


Then  for  almost  all  a  £  int  1B;iM  satisfying  Ui  —  oq  >  K2(iq  —  If)/ 2  for  1  <  i  <  n  and 
pa  defined  as  in  Theorem  4.1,  there  is  a  zero  curve  ya  of  pa  emanating  from  (0,a), 
along  which  V/Ja( A,  x)  has  full  rank,  that  remains  in  [0, 1)  x  (int  and  reaches  a 

point  (1,5:),  where  x  solves  the  MCP.  ya  does  not  intersect  itself,  is  disjoint  from  any 
other  zeros  of  pa,  and  has  finite  arc  length  if  F  is  strongly  regular  at  x. 

Proof.  The  first  four  hypotheses  of  Proposition  2.2  are  satisfied  trivially.  The 
choice  of  p(X),  and  the  restrictions  on  a  suffice  to  carry  out  the  proof  of  Theo¬ 
rem  4.1.  Hence  ya  remains  in  [0, 1)  x  (int  and  is  bounded  since  is  bounded. 

□ 

The  remainder  of  this  section  generalizes  the  above  results  to  other  choices  of  <j> 
and  (ffj,. 

Theorem  4.3.  Let.  <p  be  a  positively  oriented  NCP  function,  and  let  4>  be  a 
C2 -smoother  for  <j>,  monotone  in  its  first  two  variables,  satisfying 

(4.6)  <t>(a,/3)  >  <p(a,P,p),  for  all  a,  f3  £\ R,  p>  0,  and 


~  QljP 

(4.7)  </>(a,  0,  p)  > - ,  for  p  >  0,  0  <  a  <  oo, 

a 

where  c  andp  are  positive  constants.  Define  if  by  (2.8)  and  the  smoother  if  by  (2.12). 
Choose  a  £  int  IB/, u<  and  let  p  :  [0, 1]  — >  ]R+  be  a  decreasing  C 2  function  with  p(  1)  =  0 
satisfying 

(4.8)  p{ \)p  <  ^ c^iui  -  af){ui  -  If)  for  i  £  A  £  (0, 1] . 

Define  by  (2.14),  define  pa  '■  [0, 1)  xJR,"  — >  ]R"  by  (3.2),  and  let"ja  be  the  connected 
component  of  pf1^^})  containing  (0,a).  Then  ya  is  contained  in  [0,1)  x  (intlB;^)- 
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Proof.  The  proof  is  identical  to  the  proof  of  Theorem  4.1  except  that  in  place  of 
(4.4),  the  inequality  (4.7)  is  used.  Then  by  similar  arguments  using  (4.8), 


F^X\x)>- 


Cfl‘ 


p 


Ui  li 

1- A, 


>  — k 

A 


ai), 


contradicting  (4.3).  □ 

Corollary  4.4.  Let  cp,  cp,  ip,  ip,  and  F M  be  defined  as  in  Theorem  4.3.  Assume 
that  all  the  bounds  of  the  MCP  are  finite,  choose  k  £  (0,1),  and  take 


H(\)  =  K 


1  —  A 


i  /p  r 


min(rtj 

.  i 


Then  for  almost  all  a  £  int  satisfying  Ui  —  a,  >  kp(u,  —  U)  for  1  <  i  <  n  and  pa 
defined  as  in  Theorem  4.1,  there  is  a  zero  curve  y0  of  pa  emanating  from  (0,a),  along 
which  V pa (A,  x)  has  full  rank,  that  remains  in  [0,1)  x  (int  1B;jU)  and  reaches  a  point 
(1,5),  where  x  solves  the  MCP.  ja  does  not  intersect  itself  is  disjoint  from  any  other 
zeros  of  pa,  and  has  finite  arc  length  if  F  is  strongly  regular  at  x. 

4.1.  Ensuring  feasibility.  Since  some  MCP  applications  involve  functions  that 
are  not  defined  outside  the  feasible  region,  the  algorithm  includes  an  option  to  ensure 
that  all  iterates  are  feasible.  The  following  discussion  assumes  that  the  MCP  algorithm 
is  based  on  the  particular  choices  of  (p  and  (p^  given  by  (2.6)  and  (2.11). 

Feasibility  of  the  path  7a  can  be  assured  by  Theorem  4.1,  provided  that  the  initial 
point  a  and  the  function  p{\)  are  chosen  appropriately.  The  following  procedure 
achieves  this  while  choosing  the  initial  point  a  near  the  starting  point  x°  provided  by 
the  user:  define  a  by 


(4.9)  ai 


mid^i  +  Vi,x\,  Ui  -  vf), 
ma x(li  +  v,  x?) 
min  (iii  —  v,x®) 


if  i  €  Ipu, 

for  Ui  =  oo,  h  finite, 
for  li  =  —oo ,Ui  finite, 
if  li  =  —oo,  Ui  =  oo, 


where  17  :=  ~~  ^)/2  f°r  i  G  h,u,  and  Kmin  G  (0,1)  and  v  >  0  are  constants 

that  ensure  the  strict  feasibility  of  a.  Next,  define  p( A)  by  (3.1),  with  a  given  by 


(4.10) 


(  min  (c,  k  [minie7i  u  (ut  -  If)] )  ,  if  Ipu  ^  0, 
1  c  otherwise, 


where  c  is  some  positive  constant,  and  k  is  defined  as  follows  if  Ii,v  ^  0: 


(4.11) 


k  :=  min 

ieh.u 


Note  that  if  Ipu  is  not  empty,  this  choice  of  a  and  k  ensures  that  k  >  Kmin  and 
also  that  (4.1)  is  satisfied.  Thus,  the  assumptions  of  Theorem  4.1  are  satisfied,  so 
70  remains  strictly  within  the  feasible  region.  Feasibility  is  maintained  by  exploit¬ 
ing  STEPNX’s  built-in  logic  for  handling  domain  violations.  Precisely,  whenever  a 
STEPNX  call  to  evaluate  F(x)  produces  an  infeasible  point  (either  in  the  prediction 
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phase  or  the  correction  phase),  that  domain  violation  is  reported  to  STEPNX.  The 
result  is  that  STEPNX  cuts  the  step  size  in  half  (after  sanity  checks  to  prevent  an 
infinite  loop)  and  calculates  a  new  predicted  point.  Since  the  zero  curve  is  strictly 
feasible  for  A  <  1,  eventually  (assuming  adequate  machine  precision)  a  feasible  step 
will  be  taken. 

Finally,  to  ensure  feasibility  of  the  iterates  generated  in  the  end  game,  the  gen¬ 
eralized  damped  Newton  method  in  Figure  2.1  is  modified  according  to  the  general 
descent  framework  described  in  [12].  Specifically,  the  Newton  direction  dk  is  projected 
back  onto  the  feasible  region  to  produce  the  modified  direction 

dk  ■=  nmi  u(xk  +  dk)  -xk. 


Note  that  xk  +  dk  is  feasible.  Step  3  in  Figure  2.1  is  then  replaced  with  the  following: 
Step  3’  If  6{xk  +  dk)  <  (1  —  cr)9(xk),  set  xk+1  =  xk  +  dk.  Otherwise,  take  a 
projected  gradient  step  as  follows:  let  mk  be  the  smallest  nonnegative 
integer  ?n  <  mmax  such  that 

(4.12)  9{xk{am))  <  9(xk)  -  <7V0(: xk)(xk  -  xk(am)), 


where  xk(t)  :=  u(xk  —  tS79(xk)).  If  no  such  mk  exists,  stop;  the 
algorithm  failed.  Otherwise,  set  xk+1  =  xk(ak). 


Note  that  for  any  feasible  x* 


+  dk  -  x* 


<  \\xk  +  dk  —  . 


This  ensures. 


by  [12,  Theorem  4.5]  and  Theorem  2.1,  that  in  a  neighborhood  of  a  strongly  regular 
solution  x,  the  iterates  generated  by  the  feasible  end  game  strategy  described  above 
converge  Q-superlinearly  to  x. 

The  projected  gradient  step  in  the  above  algorithm  requires  that  9  be  differen¬ 
tiable.  This  is  true  when  <fr  is  the  Fischer-Burmeister  function  (2.6),  but  is  not  true 
in  general. 


5.  Solver  implementation  and  testing.  The  MCP  algorithm  described  in  the 
previous  section  was  implemented  using  the  Fischer-Burmeister  NCP  function  for  <f> 
and  the  smoother  defined  by  (2.11).  The  nonsmooth  Newton’s  method  described  in 
Figure  2.1  was  used  for  the  Newton  end  game.  To  construct  the  homotopy  mapping 
defined  in  (3.2),  the  parameter  a  was  constructed  according  to  (4.9),  with  Kmjn  :=  0.1, 
and  v  =  0.0001.  The  function  fi( A)  was  defined  by  (4.10),  with  c  =  1.0,  and  k  defined 
by  (4.11). 

The  algorithm  was  implemented  in  C  with  a  link  to  the  Fortran  90  subroutine 
STEPNX  from  HOMPACK90.  The  code  is  interfaced  with  the  GAMS  modeling 
language,  enabling  it  to  be  tested  using  the  MCPLIB  suite  of  GAMS  test  problems 
[11,  4].  All  linear  algebra  was  performed  using  the  LUSOL  sparse  factorization  routine 
[14]  from  MINOS  [17]. 

Computational  results  on  the  MCPLIB  problems  are  shown  in  Table  5.1.  Many 
of  the  problems  in  this  test  library  include  multiple  runs,  which  vary  the  starting 
point  x°  or  other  parameters  defining  the  problem.  All  of  the  problems  were  run 
using  default  parameter  settings,  and  the  number  of  successes  and  failures  over  all 
runs  are  reported  in  the  third  column  of  Table  5.1.  The  notation  m(n)  means  that 
the  problem  included  m  +  n  runs,  and  for  those,  there  were  m  successes  and  n  failures. 
The  default  parameters  were  chosen  as  follows: 

•  Curve  tracking  parameters:  abserr  =  relerr  =  ICC4.  Maximum  step  size 
ftmax  =  100,000.  The  normal  default  for  this  parameter  used  by  HOM- 
PACK90  is  Vax  =  1-  However,  many  problems  in  the  MCPLIB  test  library 
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were  poorly  scaled  so  had  very  long  zero  curves.  The  large  value  of  /imax  was 
therefore  used  to  allow  these  curves  to  be  tracked  in  a  reasonable  number  of 
iterations.  All  other  curve  tracking  parameters  were  the  defaults  chosen  by 
STEPNX. 

•  Newton  parameters  (See  Figure  1):  a  =  c r  =  0.5.  ramax  =  20.  Maximum 
number  of  Newton  iterations  =  30. 

•  Stopping  criteria:  An  iterate  xk  was  considered  to  solve  the  problem  when 

lljr(a;fc)J!0o  /  (1  -  <  10_6- 

In  cases  where  the  problem  was  not  solved  by  the  default  parameters,  the  algo¬ 
rithm  was  restarted  using  more  conservative  parameters:  abserr  =  relerr  =  10-6, 
dideal  =  0.01,  lideal  =  0.01,  rideal  =  0.005,  and  /imax  =  max(.l,  arclen/100), 
where  arclen  is  the  arc  length  of  the  zero  curve  calculated  using  the  default  param¬ 
eters.  Results  from  these  runs  are  shown  in  the  fourth  column  of  Table  5.1. 

For  the  problems  that  were  not  solved  by  the  conservative  settings,  the  last  column 
of  Table  5.1  describes  the  reason  for  failure.  The  notation  “oo”  indicates  that  the 
zero  curve  appeared  to  go  off  to  infinity.  This  behavior  is  common  for  problems  that 
do  not  satisfy  the  global  monotonicity  assumption.  The  notation  “lost”  indicates 
that  STEPNX  was  unable  to  continue  tracking  the  zero  curve.  This  is  generally 
due  to  a  poorly  conditioned  Jacobian  matrix.  The  notation  “r”  indicates  failure  due 
to  exceeding  resource  limits-either  the  limit  of  5000  homotopy  steps,  or  1000  CPU 
seconds.  Finally,  the  notation  “v”  indicates  failure  due  to  domain  violations. 

While  the  algorithm  failed  to  solve  a  number  of  problems  that  have  been  solved 
by  other  algorithms,  it  is  encouraging  to  note  that  it  performed  very  well  on  some 
problems  that  are  generally  regarded  as  very  hard.  Notable  among  these  are  the 
billups,  pgvonl05,  pgvonl06,  and  simple-ex  problems.  Thus,  the  homotopy  algorithm 
should  be  viewed  as  an  important  supplement  to  other  approaches. 

It  should  also  be  noted  that  the  algorithm  solved  several  problems  for  which  it  was 
not  able  to  track  the  zero  curve  all  the  way  to  A  =  1.  This  occurred  for  the  bert_oc, 
obstacle,  opt_cont*  problems.  However,  for  these  problems  the  Newton  end-game  was 
able  to  find  the  solution. 

Except  for  the  cases  “v”  and  “r”,  the  failures  are  of  two  types:  numerical  instabil¬ 
ity  or  unbounded  homotopy  zero  curve  ja.  No  attempt  was  made  to  scale,  reformulate, 
or  precondition  the  test  problems,  or  to  tune  the  tracking  parameters  for  a  particular 
problem.  There  is  little  doubt  that  a  concerted  pursuit  of  all  these  options  would 
have  removed  all  the  failures  due  to  numerical  instability.  The  unbounded  zero  curves 
are  a  more  fundamental  problem,  indicating  that  the  default  homotopy  map  (3.2) 
is  inadequate  (which  is  no  surprise,  since  in  engineering  practice  the  default  map  is 
virtually  never  used).  It  is  likely  that  replacing  (3.2)  by  XF^x>(x)  +  (1  —  A )G(a,  A,  x), 
where  G  is  carefully  crafted  for  each  problem,  could  remove  the  other  failures.  This 
remains  the  topic  of  future  work. 

6.  Conclusions.  This  paper  described  a  probability-one  homotopy  algorithm 
for  solving  nonsmooth  systems  of  equations  and  complementarity  problems.  These 
methods  are  an  extension  to  nonsmooth  equations  of  the  probability-one  homotopy 
methods  described  in  [8,  21,  23,  24]  and  they  are  attractive  because  they  are  able 
to  solve  a  qualitatively  different  class  of  problems  than  methods  relying  on  merit 
functions.  This  claim  is  justified  both  theoretically  and  computationally.  The  key 
to  success  of  the  method  is  the  global  monotonicity  assumption.  When  this  is  sat¬ 
isfied,  the  zero  curve  is  known  to  lead  to  a  solution.  This  result  is  formalized  in 
Theorem  3.2.  In  the  case  of  complementarity  problems,  an  easily  satisfiable  condition 
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Table  5.1 

MCPLIB  Test  Problems 


Problem 

Name 

size 

Default  Settings 
Success(Failure) 

Conservative  Settings 
Success  (Failure) 

Notes 

badfree 

5 

1(0) 

bert.oc 

5000 

3(1) 

3(1) 

r 

bertsekas 

15 

5(1) 

6(0) 

billups 

1 

3(0) 

bratu 

5625 

1(0) 

choi 

13 

1(0) 

colvdual 

20 

4(0) 

colvnlp 

15 

6(0) 

colvtemp 

20 

4(0) 

cycle 

1 

1(0) 

degen 

2 

1(0) 

duopoly 

63 

0(1) 

0(1) 

oo 

ehl_k40 

41 

2(1) 

3(0) 

ehl_k60 

61 

2(1) 

3(0) 

ehl_k80 

81 

2(1) 

3(0) 

ehLkost 

101 

1(2) 

1(2) 

lost 

electric 

158 

0(1) 

0(1) 

oo 

eta2100 

296 

0(1) 

1(0) 

explcp 

16 

1(0) 

forcebsm 

184 

0(1) 

0(1) 

oo 

forcedsa 

186 

0(1) 

0(1) 

oo 

freebert 

15 

7(0) 

gafni 

5 

3(0) 

games 

16 

25(0) 

hanskoop 

14 

10(0) 

hydroc06 

29 

0(1) 

0(1) 

oo 

hydroc20 

99 

0(1) 

0(1) 

oo 

jel 

6 

2(0) 

josephy 

4 

8(0) 

kojshin 

4 

8(0) 

lincont 

419 

0(1) 

0(1) 

oo 

mathinum 

3 

6(0) 

mathisum 

4 

7(0) 

methan08 

31 

0(1) 

0(1) 

oo 

multi-v 

48 

0(3) 

0(3) 

lost 

nash 

10 

4(0) 

ne-hard 

3 

1(0) 

obstacle 

2500 

7(1) 

8(0) 

was  established,  which  ensures  that  the  homotopy  zero  curve  always  remains  strictly 
feasible.  This  condition  can  always  be  enforced  in  the  algorithm  by  choosing  the  ini¬ 
tial  point  a  properly.  A  simple  consequence  of  this  result  is  that  for  finitely  bounded 
mixed  complementarity  problems,  the  zero  curve  is  bounded,  and  by  Proposition  2.2, 
is  guaranteed  to  lead  to  a  solution. 
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Table  5.1 

MCPLIB  Test  Problems  (coni.) 


Problem 

Name 

size 

Default  Settings 
Success(Failure) 

Conservative  Settings 
Success(Failure) 

Notes 

olg 

249 

0(1) 

0(1) 

lost 

opt_contl27 

4096 

1(0) 

opt_cont 

288 

1(0) 

opt_cont255 

8192 

1(0) 

opt_cont31 

1024 

1(0) 

opt_cont511 

16384 

1(0) 

pgvonl05 

105 

4(0) 

pgvonl06 

106 

5(1) 

6(0) 

pies 

42 

0(1) 

1(0) 

powell 

16 

5(1) 

5(1) 

oo 

powelLmcp 

8 

6(0) 

qp 

4 

1(0) 

romer 

214 

0(2) 

0(2) 

lost 

scarbsum 

40 

1(1) 

2(0) 

scarfanum 

13 

4(0) 

scarfasum 

14 

1(3) 

1(3) 

V 

scarfbnum 

39 

0(2) 

2(0) 

scarfbsum 

40 

1(1) 

2(0) 

shubik 

30 

7(41) 

13(35) 

r 

simple-ex 

17 

1(0) 

simple-red 

13 

1(0) 

sppe 

27 

3(0) 

tinloi 

146 

10(54) 

64(0) 

tobin 

42 

4(0) 

trade  12 

600 

1(1) 

1(1) 

lost 

trafelas 

2376 

0(2) 

0(2) 

r 
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